#include "hw1.h"
#include <algorithm>
#include <random>
#include <iostream>
using namespace algebra;

Matrix algebra::zeros(std::size_t n, std::size_t m) {
    if( n == 0 || m == 0) {
        return  Matrix{};
    }
    return Matrix(n,std::vector<double>(m,0));
}

Matrix algebra::ones(std::size_t n, std::size_t m) {
    if( n == 0 || m == 0) {
        return  Matrix{};
    }
    return Matrix(n,std::vector<double>(m,1.f));
}

Matrix algebra::random(std::size_t n,std::size_t m,
                        double min, double max) {
    if( n == 0 || m == 0) {
        return  Matrix{};
    }
    if( max < min ) {
        throw std::logic_error("In random : min cannot be bigger than max");
    }
    std::random_device rd;
    std::default_random_engine gen(rd());
    std::uniform_real_distribution<double> dis(min,max);

    auto dgen = [&]{ double temp = dis(gen); return temp; };
    
    Matrix randmtx = Matrix(n,std::vector<double>(m));
    for( auto &v : randmtx ) {
        std::generate(v.begin(), v.end(), dgen);
    }
    
    return randmtx;
}

void algebra::show(const Matrix& mat) {
    if(mat.empty()) {
        std::cout << "<empty>" << std::endl;
        return;
    }
    else if(mat[0].empty()) {
        std::cout << "<empty>" << std::endl;
        return;
    }
    for(const auto& row : mat) {
        for(const auto& ele : row) {
            std::cout << " " << ele;
        }
        std::cout << std::endl;
    }
}

Matrix algebra::multiply(const Matrix &matrix, double c) {
    if( matrix.empty() )
        return Matrix{};
    if( matrix[0].empty())
        return Matrix{};
    
    std::size_t m,n;
    m = matrix.size();
    n = matrix[0].size();

    if(c == 0.f ) {
        return Matrix(m,std::vector<double>(n,0.f));
    }
    Matrix newone(matrix.begin(),matrix.end());
    for(auto &v : newone) {
        for( auto &ele : v) {
            ele *= c;
        }
    }
    return newone;
};

Matrix algebra::multiply(const Matrix& matrix1, const Matrix& matrix2) {
    // case 1 : both empty
    if( matrix1.empty() && matrix2.empty()  ) {
        return Matrix{};
    }
    // case 2 : only one is empty
    else if(matrix1.empty() || matrix2.empty()) {
        throw std::logic_error("Not sure");
    }
    else {
        std::size_t m1 , n1 ;
        std::size_t m2 , n2 ;

        m1 = matrix1.size();n1 = matrix1[0].size(); 
        m2 = matrix2.size();n2 = matrix2[0].size();
        if( n1 != m2 || n1 == 0 || n2 == 0) {
            throw std::logic_error("deminsion error");
        }
        
        Matrix factor(m1,std::vector<double>(n2));

        for(std::size_t i = 0; i < m1 ; i ++) {
            const auto& r_now = matrix1[i];
            
            for(std::size_t j = 0; j < n2 ; j++) {
                double vec_factor = 0;
                for(std::size_t k = 0; k < n1 ; k++) {
                    vec_factor += r_now[k] * matrix2[k][j];
                }
                factor[i][j] = vec_factor;
            }
        }
        return factor;
    }
}

Matrix algebra::sum(const Matrix &matrix, double c) {
    if(matrix.empty()) {
        return Matrix{};
    }

    std::size_t m,n;
    m = matrix.size();
    if( matrix[0].empty() )
        return Matrix{};
    n = matrix[0].size();
    
    Matrix newone = matrix;

    for(auto & v : newone) {
        for (auto & ele : v) {
            ele += c;
        }
    }

    return newone;
}

Matrix algebra::sum(const Matrix& matrix1, const Matrix& matrix2) {
    if(matrix1.empty() && matrix2.empty()) {
        return  Matrix{};
    }
    else if(matrix1.empty() || matrix2.empty()) {
        throw std::logic_error("Not sure");
    }
    else {
        std::size_t m1,n1;
        std::size_t m2,n2;

        m1 = matrix1.size(); n1 = matrix1[0].size();
        m2 = matrix2.size(); n2 = matrix2[0].size();

        if(m1 != m2 || n1 != n2 ) {
            throw std::logic_error("dimension not match");
        }
        Matrix ret(m1,std::vector<double>(n1));
        for(std::size_t i = 0 ; i < m1 ; i ++) {
            for(std::size_t j = 0; j < n1 ; j ++) {
                ret[i][j] = matrix1[i][j] + matrix2[i][j];
            }
        }
        return ret;
    }
}

Matrix algebra::transpose(const Matrix &matrix) {
    if(matrix.empty()) {
        return Matrix{};
    }
    else {
        std::size_t m = matrix.size();
        std::size_t n = matrix[0].size();

        if( n == 0) {
            throw std::logic_error("Not Sure");
        }

        Matrix ret(n,std::vector<double>(m));
        for(std::size_t i = 0 ; i < n ; i++) {
            for(std::size_t j = 0 ; j < m ; j ++) {
                ret[i][j] = matrix[j][i];
            }
        }
        return ret;


    }
}
Matrix algebra::minor(const Matrix &matrix, std::size_t n, std::size_t m) {
    if( matrix.empty() ) {
        throw std::logic_error("Matrix is empty");
    }
    else {
        std::size_t m_row, m_col;
        m_row = matrix.size();
        m_col = matrix[0].size();

        //should not be empty
        if(m_col == 0) {
            throw std::logic_error("Matrix is empty");
        }

        //should be a square and not a 1x1 square
        if(m_col != m_row ) {
            throw std::logic_error("Matrix is not a square");
        }
        if(m_col == 1 || m_row == 1) {
            throw std::logic_error("Matrix is a 1x1 square");
        }
        // nth row , m th col will be removed
        Matrix ret(m_row - 1, std::vector<double>(m_col - 1));
        std::size_t  minor_i, m_i;
        std::size_t  minor_j, m_j;
        for( minor_i = 0,m_i = 0; minor_i < (m_row - 1);) {
            if(m_i != n) {
                for(minor_j = 0,m_j = 0 ; minor_j < (m_col - 1);) {
                    if(m_j != m) {
                        ret[minor_i][minor_j] = matrix[m_i][m_j];
                        minor_j++;
                    }
                    m_j ++;
                }
                minor_i ++;
            } 
            m_i ++;
        }
        return ret;
    }
}

double algebra::determinant(const Matrix &matrix) {
    if(matrix.empty()) {
        // throw std::logic_error("Matrix is empty");
        return 1;
    }
    else {
        std::size_t m,n;
        m = matrix.size();
        n = matrix[0].size();

        //should not be empty and should be square
        if(n == 0)
            throw std::logic_error("Matrix is empty");
        if(n != m)
            throw std::logic_error("Matrix is not a square");

        // now sure it's a square
        if(m == 1 && n == 1) 
            return matrix[0][0];
        else {
            // choose 1st row to calculate
            int factor = 1;
            double accu = 0.0f;
            for(std::size_t j = 0; j < n; j++) {
                // a_i_j * C_i_j
                double a = matrix[0][j];
                double C = factor * algebra::determinant( algebra::minor(matrix,0,j));

                accu += a * C;
                factor = -factor;
            }
            return accu;
        }
    }
}
Matrix algebra::inverse(const Matrix& matrix) {
    if(matrix.empty()) {
        return Matrix{};
    }
    else {
        std::size_t m,n;
        m = matrix.size();
        n = matrix[0].size();

        //should not be empty and should be square
        if(n == 0)
            throw std::logic_error("Matrix is empty");
        if(n != m)
            throw std::logic_error("Matrix is not a square");

        // now sure it's a square
        if( m == 1) {
            double elem = matrix[0][0];
            if( std::fabs(elem - 0.f) < 1e-6 )
                throw std::logic_error("singular Matrix has no inverse");
            else 
                return Matrix(1,std::vector<double>(1, 1.f / (elem)));
        }
        else {
            // Cal det -> Cal adj ;
            // Step 1 : calculate adj

            double det = determinant(matrix);
            if( std::fabs(det - 0.f)  < 1e-6 )
                throw std::logic_error("singular Matrix has no inverse");
            
            // Step 2 : calculate adj , then directly get the reverse
            std::size_t m,n;
            m = matrix.size();
            n = matrix[0].size();

            Matrix rev(m,std::vector<double>(n));
            int factor = 1;
            for(std::size_t i = 0 ; i < m ; i ++ ) {
                for(std::size_t j = 0; j < n ; j++) {
                    double c = algebra::determinant( algebra::minor(matrix, i , j) );
                    c = factor * c;
                    rev[j][i] = c / det;
                }
            }
            return rev;
        }
    }
}

Matrix algebra::concatenate(const Matrix &matrix1, const Matrix &matrix2 ,int axis = 0) {
    if(matrix1.empty() || matrix2.empty() ) {
        throw std::logic_error("One of Matriecs is empty");
    }
    if( axis == 0 ) {
        std::size_t m1,n1;
        std::size_t m2,n2;

        m1 = matrix1.size(); n1 = matrix1[0].size();
        m2 = matrix2.size(); n2 = matrix2[0].size();

        // should not be empty and demision matches
        if( n2 == 0 || n1 == 0 ) {
            throw std::logic_error("One of Matriecs is empty");
        }
        if( n2 != n1) {
            throw std::logic_error("demision not match");
        }

        Matrix ret(m1 + m2 , std::vector<double>(n1));
        for(std::size_t i = 0 ; i < (m1 + m2) ; i++) {
            for(std::size_t j = 0; j < n1 ; j ++) {
                ret[i][j] = i >= m1 ? matrix2[i - m1][j] : matrix1[i][j];
            }
        }
        return  ret;
    }
    else {
        std::size_t m1,n1;
        std::size_t m2,n2;

        m1 = matrix1.size(); n1 = matrix1[0].size();
        m2 = matrix2.size(); n2 = matrix2[0].size();

        // should not be empty and demision matches
        if( n2 == 0 || n1 == 0 ) {
            throw std::logic_error("One of Matriecs is empty");
        }
        if( m1 != m2) {
            throw std::logic_error("demision not match");
        }

        Matrix ret(m1, std::vector<double>(n1 + n2));
        for(std::size_t i = 0 ; i < (m1) ; i++) {
            for(std::size_t j = 0; j < (n1 + n2) ; j ++) {
                ret[i][j] = j >= n1 ? matrix2[i][j - n1] : matrix1[i][j];
            }
        }
        return ret;
    }
}

Matrix algebra::ero_swap(const Matrix &matrix, size_t r1, size_t r2) {
    if( matrix.empty() ) 
        throw std::logic_error("Matrix is empty");
    else {
        std::size_t m,n;
        m = matrix.size();
        n = matrix[0].size();

        if(n == 0) 
            throw std::logic_error("Matrix is empty");

        if( r1 >= m || r2 >=m )
            throw std::logic_error("row number out of range");

        Matrix ret{matrix};
        std::swap(ret[r1],ret[r2]);

        return ret;
    }
}

Matrix algebra::ero_multiply(const Matrix &matrix, size_t r, double c) {
    if( matrix.empty() )
        throw std::logic_error("Matrix is empty");
    else {
        std::size_t m,n;
        m = matrix.size();
        n = matrix[0].size();

        if( n == 0) 
            throw std::logic_error("Matrix is empty");

        if( r >= m )
            throw std::logic_error("row number out of range");

        if( fabs(c)  < 1e-10)
            throw std::logic_error("ERO multiply not allow zero multiply");
        Matrix ret{matrix};
        for(double &elem : ret[r]) {
            elem *= c;
        }

        return ret;
    }
}

Matrix algebra::ero_sum(const Matrix &matrix, size_t r1, double c, size_t r2) {
    if( matrix.empty() )
        throw std::logic_error("Matrix is empty");
    else {
        std::size_t m,n;
        m = matrix.size();
        n = matrix[0].size();

        if(n == 0) 
            throw std::logic_error("Matrix is empty");

        if(r1 >= m || r2 >=m)
            throw std::logic_error("row number out of range");

        if(fabs(c) < 1e-10)
            throw std::logic_error("ERO multiply not allow zero multiply");
        
        
        Matrix ret{matrix};
        for(std::size_t i = 0; i < n ; i ++)
            ret[r2][i] += ret[r1][i] * c;

        return ret;
    }
}

Matrix algebra::upper_triangular(const Matrix &matrix) {
    if( matrix.empty() )
        return Matrix{};
    else {
        std::size_t m,n;
        m = matrix.size();
        n = matrix[0].size();

        // should not be empty one
        if(n == 0)
            return Matrix{};
        // should be a square
        if( m != n )
            throw std::logic_error("matrix is not a square.");

        // now it's a square
        if( m == 1 ) {
            return  Matrix{matrix};
        }
        else {
            Matrix ret{matrix};
            std::size_t major = 0;
            std::size_t row   = 0;
            for(major = 0; major < m ;) {
                // check the position is bigger than 0
                if( fabs(ret[row][major]) <= 1e-10 ) {
                    // if it's 0, consider to find another row.
                    std::size_t non_zero_row;
                    for( non_zero_row = row + 1 ; 
                         non_zero_row < m; 
                         non_zero_row++
                    ) {
                        if( fabs(ret[non_zero_row][major]) > 1e-10) {
                            break;
                        }
                    }
                    
                    if( non_zero_row < m ) {
                        // find and switch with the row now
                        std::swap(ret[non_zero_row],ret[row]);
                    }
                    else {
                        // all zero here , find next major
                        major ++;
                        continue;
                    }

                }
                // for rest lines , perform ERO mul-add
                for( size_t i = row + 1 ; i < m ; i++ ) {
                    // check the major of this row is zero
                    if( fabs(ret[i][major]) > 1e-10  ) {
                        // if not , Perform ERO add-mul operation
                        double factor = -(ret[i][major] / ret[row][major]) ;
                        for(size_t j = major + 1 ; j < m ; j++) {
                            ret[i][j] += factor * (ret[row][j]);
                        }
                    }
                    ret[i][major] = 0;
                }
                // find next major , begin with next line
                major ++;
                row ++;
            }
            return ret;
        }
    }

}